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Nanoscale physics and dynamical mean field theory have both generated increased interest in complex quan- 
tum impurity problems and so have focused attention on the need for flexible quantum impurity solvers. Here 
we demonstrate that the mapping of single quantum impurity problems onto spin-chains can be exploited to 
yield a powerful and extremely flexible impurity solver. We implement this cluster algorithm explicitly for the 
Anderson and Kondo Hamiltonians, and illustrate its use in the "mesoscopic Kondo problem". To study univer- 
sal Kondo physics, a large ratio between the effective bandwidth D c h and the temperature T is required; our 
cluster algorithm treats the mesoscopic fluctuations exactly while being able to approach the large D c ff /T limit 
with ease. We emphasize that the flexibility of our method allows it to tackle a wide variety of quantum impurity 
problems; thus, it may also be relevant to the dynamical mean field theory of lattice problems. 
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The Kondo problem, which describes the coupling of a 
quantum magnetic impurity to an electron gas, has a rich his- 
tory in condensed matter physics. 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 Recent experi- 
mental and theoretical developments have generated consid- 
erable renewed interest in this field. On the experimental 
side, nanostructures in ultra-small metal particles, semicon- 
ductor heterostructures, carbon nanotubes, and organometal- 
lic molecules all show physics that can be mapped onto quan- 
tum impurity p ro blemsi2ii2iIiiI2ii2ii2iI2iil!iii On the theoretical 
side, dynamical mean field theory (DMFT) has reduced lattice 
problems in infinite dimensions to self-consistent quantum 
impurity problems. 18 Nanoscale physics and DMFT have fur- 
thermore emphasized the richness and variety of these ques- 
tions: In mesoscopic physics quantum dots behave as impu- 
rities; their tunability, as well as the new energy scales im- 
plied by the confined geometry, have lead to several new im- 
purity problems. In DMFT, lattice models map onto complex 
self-consistent quantum impurity problems that may not have 
arisen in previous contexts. Clearly, the paradigm of a spin- 
1/2 magnetic impurity in a metallic host proposed by Kondc 2 
is only a starting point in exploring a vast array of physically 
relevant quantum impurity problems. 

In order to understand quantum impurity problems in depth 
and to compare with experiments, it is desirable to have ex- 
act results. Since exact analytic solutions 4 ^ make approxi- 
mations that are invalid beyond the simplest situations, it is 
necessary to resort to numerical approaches. The need for 
precise numerical solvers has been especially highlighted by 
DMFT. 18 Many numerical methods have been developed for 
quantum impurity problems, including the numerical renor- 
malization group (NRG) 3 and the quantum Monte Carlo algo- 
rithm introduced by Hirsch and Fye. 6 Although these meth- 
ods are powerful, they have individual weaknesses related to 
the principles on which they are formulated. Here we pro- 
pose a new quantum Monte Carlo method for quantum impu- 
rity problems. Being formulated on a different principle from, 
for instance, either NRG or the Hirsch and Fye Monte Carlo 
method, this new approach has a different set of strengths and 
weaknesses; it is, in fact, more efficient than the traditional 
methods in some cases where these latter have difficulties. 

One example where the weaknesses of the standard meth- 



ods show up is the "mesoscopic Kondo problem". Nanoscale 
experiments have motivated new questions that are usually ig- 
nored in the study of quantum impurity problems. One may 
ask, for instance, what is the effect of confining the bath elec- 
trons to a finite, fully coherent, region of space. In this case, 
the local electronic density of states fluctuates strongly, and a 
flat band description is no longer applicable. In addition to the 
effects caused by a finite mean level spacing AjI2i22i2i confine- 
ment leads to mesoscopic fluctuations associated with inter- 
ference effects^S 3 . that cannot be captured exactly by analytic 
methods. Numerically, NRG is not natural for a situation with 
discrete levels, and involves some approximations (essentially 
neglecting part of the conduction electrons' Hilbert space) that 
have yet to be tested for structured density of states. Although 
the Hirsch and Fye algorithm does not depend on the con- 
duction band, the formulation of the algorithm makes it im- 
practical to use in the most interesting regime when Z3 c ff (the 
effective bandwidth where the impurity is magnetic) becomes 
large compared to the temperature T. Hence, a detailed study 
of the effects of mesoscopic fluctuations on universal Kondo 
physics requires new numerical tools, able to handle any kind 
of discrete density of states while allowing the study of large 
D c g/T. The goal of this paper is to show that such a QMC 
algorithm indeed exists. 

We develop our QMC method explicitly for two quan- 
tum impurity problems, the Anderson and Kondo models, 
and point out later that our algorithm is very flexible and 
has a broad range of applications. The Hamiltonian for 
both models can be written as a sum of two terms H = 
Hq + Hi where Ho describes the free electron band, Ho = 
2 Q =i S<t±i e Q c icr c QO'' ar, d i s entirely specified by the ener- 
gies e a and the eigenfunctions (j> a (r) of the confined electron 
gas. We will assume —D<e a < D, where D is the band- 
width. A "flat" band would correspond to the choice e a = a A 
and |0 Q (r)| 2 = I/O, which, for temperatures larger than the 
mean level spacing A would just reproduce the bulk behavior. 
Here O is the volume of the electron gas. For fully coherent 
confined systems, however, both of these quantities display 
mesoscopic fluctuations from level to level, or as a function of 
some external parameter. For chaotic systems with time rever- 
sal symmetry, a good model for these fluctuations is provided 



by the Gaussian Orthogonal Ensemble (GOE) of random ma- 
trix theory. 24 The averaged quantities are kept the same as in 
the flat band case, and in particular (p(e)) — po = l/(Af2), 
where (p(e)) is the average over realizations of the local den- 
sity of states p{e) —^2 a |^a(r)| 2 5(e — e a ). In the following, 
we shall study both a flat band and some realizations drawn 
from GOE for illustration. 

The interactions of the conduction electrons with the impu- 
rity are encoded in Hi along with the impurity Hamiltonian. 
For the Anderson modei24 

Hf = y^[*t(r K + 4vi/ ff ( ro )] (i) 

E7 

+ ^e d dld a + Ud\di;d\di. 

a 

where \fft (r ) = J2 a <$>* a {ro)c' a(T creates an electron of spin 
a at the location ro of the impurity. Similarly for the Kondo 
model 3 . 

H * = \ E ^(ro)ow^Mro) • , (2) 

where S d is the impurity spin and a is the vector formed by 
the Pauli matrices. For a flat band and when e& = —U/2 
(symmetric case) the Anderson model turns into the Kondo 
model when U — > oo with Jp — 8T/irU and D fixed. Here 
we have defined T = irpoV 2 . 

The basic strategy we are going to apply here relies on the 
fact that the single-impurity Hamiltonians are essentially one- 
dimensional problems. Indeed, for both the Anderson and 
Kondo models, the impurity couples to the electron gas only 
locally, at ro- Starting from ji„ = ^ (T (ro), one can find a 
basis of one particle states fi a ,i — l, 2..,N such that the non- 
interacting Hamiltonian Hq becomes tridiagonal 1 - 3 ' 26 - 27 : 

N 

(7 i—1 

with /3tv = 0. In this form we see that Hq describes two 
open fermion chains. Each open fermion-chain is identi- 
cal to a spin-chain since fermions hopping in one dimension 
cannot permute. Thus, there is no fermion sign problem to 
worry about when constructing quantum Monte Carlo algo- 
rithms. Adding the interaction term Hi merely couples the 
two fermion-chains at i = 1 with the impurity (see, for ex- 
ample, the illustration in Ref. 113). In this form, one merely 
has to solve a spin-chain problem. Today spin-chain prob- 
lems can be solved very efficiently in continuous time 2 ^ 2 ^ us- 
ing the recently developed directed-loop cluster algorithm. 30 
In the present case, to make the algorithm efficient one must 
perform two types of directed-loop updates: One that changes 
the fermion occupation numbers and the other that flips the 
fermion spins (changes a). Since the directed-loop algorithm 
is well established, we will not discuss it here. 

We have implemented this algorithm for both the Ander- 
son and the Kondo models using the continuous-time path 
integral directed-loop algorithm. 30 Although the algorithm is 
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FIG. 1: Comparison of the local susceptibilities \ obtained using 
the spin-chain algorithm (circles) and the Hirsch and Fye algorithm 
(open squares) in the symmetric Anderson model for N = 5000, 
2D = 40, and U/F = 8/1.6 (Tk = 0.146). %Tk is in units of 
(<7/ib/2) 2 . Three different realizations are shown: (a) Clean case 
and (b) two realizations drawn from GOE. (c) The real and imagi- 
nary parts of the thermal Green function versus Matsubara frequen- 
cies for the clean case at T/Tk = 0.55. (d) The time required in the 
spin-chain algorithm to obtain two percent errors on x as a function 
Tk/T. In the low temperature regime, an almost quadratic depen- 
dence in 1/T is observed (compared to 1/T for a single sweep at 
fixed At in the Hirsch and Fye algorithm). 6 The value of Tk used 
was obtained from the two-loop RG equation of the Kondo model. 



equally applicable to the symmetric and the asymmetric An- 
derson models, here we focus only on the symmetric case for 
convenience. The two directed-loop updates discussed above 
can readily give two correlation functions: the impurity one- 
particle thermal Green function Gd(r) and the local suscep- 
tibility x- G<j(r) for the impurity can be measured while 
directed-loops for changing the fermion occupation numbers 
are being constructed. 31 Moreover, the Green function at 

Matsubara frequencies, G d (iw„) = f^dr (d(r)<& (0)}e" J » T , 
necessary in a DMFT calculation can be directly mea- 
sured since the Fourier integral involved in Gd(«o;„) can be 
performed analytically during the loop construction in the 
continuous-time path integral implementation of the directed- 
loop algorithm. Similarly, X = (S z d (T)S z d {0)), can 
be measured while directed-loops for flipping the fermion 
spins are being constructed. We have tested our algorithm 
against exact diagonalization methods on small systems. As a 
further check, we have reproduced the results obtained from 
the Hirsch and Fye algorithm for the Anderson model. In 
Fig.^ we compare results for three different realizations of 
Hq. We see that there is perfect agreement between the Hirsch 
and Fye and the spin-chain approach in the Anderson model. 
The thermal Green function for Matsubara frequencies is also 
shown for the clean case [Fig.^c)]. 

To test our algorithm at larger U, we cannot compare to 
the Hirsch and Fye algorithm as it is very difficult to make U 
large in that approach. Hence, we compare our results to those 
of a general algorithm we recently developed to study a vari- 
ety of quantum impurity problems. 32 The local susceptibility 
obtained in the spin-chain algorithm matches very well the re- 
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FIG. 2: Local susceptibilities in the symmetric Anderson model for 
various values of U at fixed 8T/wU =Jp=l/ir, N = 1000, 2D = 10 
(Ik = 0.122) for a flat band. Data for the Kondo model at the same 
J and p is also shown. The value of Ik is obtained from the two 
loop RG equation in the Kondo model. 



suit published in Ref. 32 (parameters: N = 2000, 2D = 20, 
U = 25, and U/T = 4). In Fig.|2]we show how the Ander- 
son model results approach those of the Kondo model as U 
becomes large at fixed J AS 

We believe that the effort in the spin-chain cluster algorithm 
scales as a function of the system size as, in dimensionless 
units, either (T K /T){D/A) or (|/3i |/T)(£)/A), whichever 
dominates. This shows that there is one disadvantage of the 
spin-chain approach: Our inability to deal with a continuum 
density of states, i.e. A = 0. As one approaches this limit, the 
number of sites in the chain, N, will grow, and this in turn will 
increase the effort in updating the chain. As we have demon- 
strated, without any special effort we can simulate systems 
with N = 5000 on desktop computers; this could presumably 
be increased somewhat without too much difficulty. On the 
other hand, in the Hirsch and Fye approach the effort does not 
depend on A. We can mitigate this disadvantage by noticing 
that in the spin-chain approach it is relatively straightforward 
to apply a logarithmic blocking of the energy levels very sim- 
ilar to the one introduced by Wilson 3 in his work on the NRG. 
Similar ideas have also been developed by several authors to 
study pseudo-gap problematic in single-impurity models. 

The starting point of the logarithmic blocking is to divide 
the bandwidth [-D, D] into a finite number of energy interval 
h,h- ad I-i,I-2,... where /„ = [DA~ n , DA 1 ^ 71 ] forn > 
and I n = {-DA 1+n , -DA n ] for n < 0, with a constant A > 
1. As |n| increases, the size of a logarithmic bins becomes so 
small that there may contain only a few states. When this 
number is less than some constant I, we stop the logarithm 
discretization and keep all the remaining states. The optimal 
value for I should be determined by trial and error, but in the 
following we will just take it to be one. For each interval /„, 
we define the operator 
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where M„ = (J2 €ae i n 4> 2 a ) 1/2 - N ° te that since *J( r o) = 
^2 n Mndjjg., the interaction terms in the Anderson and Kondo 
Hamiltonians depend only on a na and a\ la . Thus, the only 
approximation in the blocking scheme comes from replacing 
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FIG. 3: Local susceptibilities using logarithmic blocking, for the 
three realizations shown in Fig. the clean case (top) and with 
mesoscopic fluctuations (bottom). The values of Na are roughly 80 
at A = 1.2, 40 at A = 1.4, and 24 at A = 2. The Kondo tempera- 
ture Tk(A) is obtained using the renormalized J e ff(A) for logarith- 
mic discretization:^ T K =0.145, 0.141, and 0.132, for A=1.2, A=1.5, 



H by H = Y^ n e 



where 



and A=2.0, respectively. 



Note that as A approaches 1, Hq —> Ho; the blocking disap- 
pears and all states are taken into account exactly. Blocking 
reduces the the number of sites necessary in the spin-chain 
formulation to A^a ~ log N. For a constant mean density of 
energy levels, A^a ~ log NJ log A which even for moderate A 
represents a very significant reduction. 

In order to check the usefulness of the logarithmic blocking 
we studied three different Hamiltonians (one clean and two 
realizations drawn from GOE). For each case we calculated 
the local susceptibilities using the blocked Hamiltonians with 
A = 1.2, 1.5, and2.0. Our results are shown in Fig. [3] One can 
see that the results from the blocked Hamiltonian approach the 
result without blocking (A = 1) as A gets closer to 1; all the 
A = 1.2 data are almost statistically indistinguishable from the 
exact results. Clearly, the optimal value of A will depend on 
the parameters of the problem, the observables measured, and 
their needed accuracies. However, the reduction in the number 
of electron sites needed is dramatic even for small A. In the 
example here, one could reduce the number of sites from 5000 
to roughly 80 with A = 1 . 2, to roughly 40 with A = 1.5, and to 
roughly 24 with A = 2.0. 

The inability to treat a continuous density of state, even 
if it can be mitigated by the logarithmic blocking described 
above, will impose some limitations on the scope of applica- 
tion of our spin-chain cluster algorithm. It has, however, a 
number of strengths which we review briefly now, and com- 
pare more particularly to the well-established Hirsch and Fye 
method. 6 Some of the advantages of the new method are: 
(1) One remarkable feature of our algorithm is illustrated in 
Fig-Ed), which shows the computer time needed to obtain the 
local susceptibility within a fixed statistical error as a func- 
tion of temperature. We find that the time grows approxi- 
mately as 1/T 2 . For the Hirsch and Fye method the time 
for a QMC sweep 6 grows as 1/T 3 , and the scaling with fixed 
fractional error is likely to be worse. Thus the present algo- 
rithm should out-perform Hirsch and Fye for low tempera- 
tures. (2) There is no imaginary time discretization error in 
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the present approach. (3) A big advantage of our spin-chain 
cluster algorithm is that it is essentially unaffected by the 
value of U. In contrast, the relevant dimensionless parameter 
determining the computational effort for the Hirsch and Fye 
Anderson-model algorithm is U/T. Since U is proportional 
to the effective bandwidth D c g within which charge fluctua- 
tion can be neglected, the ratio D e g/T is limited to a rela- 
tively modest value (~ 100) in the Hirsch and Fye algorithm. 
This value is not large enough to fully study the implications 
of the mesoscopic fluctuations on the renormalization of the 
coupling constant which dominates Kondo physics. 23 (4) We 
wish to stress that although we have illustrated our spin-chain 
approach for the two simplest impurity models, the general 
method is extremely flexible and is applicable to a wide class 
of multi-band fermionic or spin quantum impurity problems. 
The flexibility is already illustrated by the fact that the imple- 
mentation for the Kondo model Eq. (0 is essentially as simple 
as for the Anderson model Eq. (0. The Kondo version 8 of the 
Hirsch and Fye algorithm is, in contrast, substantially more 
complex than its original Anderson counterpart, and appears 
in any case significantly less used. Furthermore, although in 
Hirsch and Fye, the 1 -channel problem may not suffer from 
the sign problem, in general the 2-channel problem does. 37 
In our method there is no sign problem in either case. That 
the TV-channel Kondo problem 36 does not suffer from a sign 



problem in our method can be easily proved because, first, 
fermions in one channel cannot permute, and, second, since 
channel number is conserved, permutation of fermions in dis- 
tinct channels is also forbidden. Finally, based on our expe- 
rience we believe that the spin-chain algorithm may be ex- 
tendable to more complex impurity problems, involving for 
instance more than one orbital. 

To summarize, we have proposed a QMC method that can 
be formulated for a large class of quantum impurity problems 
in continuous time without a sign problem. We have illus- 
trated this method specifically to simulate the Anderson and 
Kondo single-impurity models, and shown how it is particu- 
larly useful for simulations of these models in the context of 
mesoscopic physics. The effort in computing quantities grows 
as a function of TkN/T. In this work N of order 5000 and 
Tk/T of the order of 10 could easily be reached. If signifi- 
cantly larger values of N are required, it is possible to approx- 
imate the problem using logarithmic blocking, similar to that 
used in the NRG approach. This reduces the effective number 
of levels in the spin-chain and should allow one to probe sig- 
nificantly lower temperature. The ideas presented here can be 
extended easily to the multi-channel Kondo model, 36 and the 
algorithm should be applicable to DMFT calculations. 18 
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* Permanent address: Laboratoire de Physique Theorique et 
Modeles Statistiques (LPTMS), 91405 Orsay Cedex, France. 

1 A.C. Hewson, The Kondo Problem to Heavy Fermions (Cam- 
bridge University Press, Cambridge, 1993). 

2 J. Kondo, Prog. Theor. Phys. 32, 37 (1964). 

3 K.G. Wilson, Rev. Mod. Phys. 47, 773 (1975). 

4 N. Andrei, Phys. Rev. Lett. 45, 379 (1980). 

5 P.B. Wiegmann, JETP Lett. 31,163 (1980). 

6 J.E. Hirsch, and R.M. Fye, Phys. Rev. Lett. 56, 2521 (1986). 

7 R.M. Fye, and J.E. Hirsch, Phys. Rev. B 38, 433 (1988). 

8 R.M. Fye, and J.E. Hirsch, Phys. Rev. B 40, 4780 (1989). 

9 T.K. Ng and PA. Lee, Phys. Rev. Lett. 61, 1768 (1988). 

10 L.I. Glazman and M.E. Raikh, JETP Lett. 47, 452 (1988). 

11 D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu, D. Abusch- 
Magder, U. Meirav and MA. Kastner, Nature 391, 156 (1998). 

12 D. Goldhaber-Gordon, J. Gores, M.A. Kastner, H. Shtrikman, 
D. Mahalu, and U. Meirav, Phys. Rev. Lett. 81, 5225 (1998). 

13 S.M. Cronenwett, T.H. Oosterkamp, and L.P Kouwenhoven, Sci- 
ence 281, 540 (1998). 

14 W.G van der Weil, S. De Franceschi, T. Fujisawa, J.M. Elzerman, 
S. Tarucha, and L.P. Kouwenhoven, Science 289, 2105 (2000). 

15 J. Nygard, D.H. Cobden, and PE. Lindelof, Nature 408, 342 
(2000). 

16 J. Park, A.N. Pasupathy, J.I. Goldsmith, C. Chang, Y. Yaish, J.R. 
Petta, M. Rinkoski, LP Sethna, H.D. Abruna, PL. McEuen, and 
D.C. Ralph, Nature 417, 722 (2002). 

17 W. Liang, M.P Shores, M. Bockrath, J.R. Long, and H. Park, Na- 
ture 417, 725 (2002). 

18 A. Georges, G. Kotliar, W. Krauth, and M.J. Rozenberg, Rev. 
Mod. Phys. 68, 13 (1996). 

19 W.B. Thimm, J. Kroha, and J. von Delft, Phys. Rev. Lett. 82, 



2143 (1999). 

20 P. Simon and I. Affleck, Phys. Rev. Lett. 89, 206602 (2002). 

21 PS. Cornaglia and C.A. Balseiro, Phys. Rev. B 66, 1 15303 (2002). 

22 R.K. Kaul, D. Ullmo, and H.U. Baranger, Phys. Rev. B 68, 
161305(R) (2003). 

23 R.K. Kaul, D. Ullmo, S. Chandrasekharan, and H.U. Baranger, 
preprint, cond-mat/040921 1 (2004). 

24 O. Bohigas, in Chaos and Quantum Physics, edited by M.J. Gi- 
annoni, A. Voros, and J. Zinn- Justin (North-Holland, Amsterdam 
1991) pp. 87-199. 

25 PW. Anderson, Phys. Rev. 124, 41 (1961). 

26 H.R. Krishna-murthy, J.W. Wilkins, and K.G. Wilson, Phys. Rev. 
B 21, 1003 (1980). 

27 C. Raas, GS. Uhrig, and F.B. Anders, Phys. Rev. B 69, 041 102(R) 
(2004). 

28 B.B. Beard and U.J. Wiese, Phys. Rev. Lett. 77, 5130 (1996). 

29 H.G Evertz, Adv. Phys. 52, 1 (2003). 

30 OF. Syljuasen and A.W. Sandvik, Phys. Rev. E 66, 046701 
(2002). 

31 A. Dorneich and M. Troyer, Phys. Rev. E 64, 066701 (2002). 

32 J. Yoo, S. Chandrasekharan, and H.U. Baranger, preprint, 
cond-mat/0408123 (2004). 

33 J.R. Schrieffer, and PA. Wolff, Phys. Rev. 149, 491 (1966). 

34 K. Chen and C. Jayaprakash, J. Physics: Cond. Matt. 7, 491 
(1995). 

35 C. Gonzalez-Buxton, and K. Ingersent, Phys. Rev. B 57, 14254 
(1998). 

36 P. Nozieres and A. Blandin, J. Phys. (Paris) 41, 193 (1980). 

37 M. Jarrell, H. Pang, D.L. Cox, and K.H. Luk, Phys. Rev. Lett. 77, 
1612 (1996). 



